In floating point units of some CPUs or, more particularly, contemporary graphics chips, it is necessary to compute mathematical functions of certain floating point numbers. These are typically continuous and may include functions such as ‘reciprocal’, i.e. ƒ(x)=1/x, ‘reciprocal square root’, ƒ(x)=1/√{square root over (x)}, ‘logarithm base 2’, ƒ(x)=ln x/ln 2, power base 2, ƒ(x)=2x, sine and cosine. For example, such mathematical functions form part of the instruction set of the Microsoft Direct X graphics specification, and some have appeared, in various forms, as instructions in various CPUs.
With many of these functions, there is a repeating pattern of behaviour that can be exploited to simplify the implementation—cosine and sine are obvious examples, but, this behaviour can also be seen with reciprocal which ‘repeats’ for each for every interval of the form [2n,2n+1) and for ‘reciprocal square root’ which repeats on intervals of the form [2n,2+2). Such regions can usually be extracted through simple integer manipulation of the input floating point exponent, and so the difficult part of the evaluation reduces to evaluating a function on a much-reduced range of numbers. Without loss of generality, it can be assumed that the input value to such a reduced function can be mapped to a (typically fixed point) value, x′ in the range [0, 1]. The output is also assumed to be of the same form.
There are several methods of computing such reduced functions known in the art. The simplest system is that of a direct table look-up. For a given input precision (N bits) and a particular output precision (M bits) the table would contain 2N entries each of M-bits. The input value is mapped to the nearest entry in the table and the result is obtained. However, for anything except small tables, this is an expensive and generally inaccurate option.
One approach that can be used to alleviate the expense is to compress the data in the table—such a scheme was used for reciprocal calculations in U.S. Pat. No. 6,330,000. This scheme, however, is difficult to apply generally and, furthermore, does not offer large levels of compression.
An alternative approach is to use a small look-up table, which only gives moderately precise results, and then apply Newton-Raphson iteration to improve the numerical accuracy. With such a scheme, the precision, in terms of numbers of bits, typically doubles with each such iteration. Turkowski describes such a method in ‘Graphics Gems V’ [ISBN 0-12-543455-3] for computing the reciprocal square root. A look-up table with, say, 26 entries, produces an initial guess, R0, which is then improved by use of two Newton-Raphson iteration steps that evaluate:Ri+1=(3−Ri*Ri*x)*x/2
Similarly, the Texas Instrument's ‘TMS320C3x User's Guide’ gives examples for the computation of ‘reciprocal’ and ‘reciprocal square root’. It, in fact, effectively dispenses with the look-up table and starts with a very approximate initial guess and simply relies on the iteration to improve the accuracy.
One problem with methods that use iterative improvements such as Newton-Raphson is that they require the additional mathematical operations (either floating point or integer) in order to converge to the result. This may not be an issue when the iterations steps are performed in a CPU with multipliers or adders that may otherwise be idle, but including the extra circuits in a dedicated hardware unit may be expensive.
Another alternative is to use a power series to evaluate the function, say, using three or four terms. In this approach, weighted versions of the k derivatives are stored as constants and the function is evaluated as . . .ƒ(x+x′)=C0+C1x′+C2x′2+ . . . +Ckx′k  Eqn. 1or, equivalently, with fewer multiplies, as . . .ƒ(x+x′)=C0+(C1+(C2+ . . . )x′)x′  Eqn 1a
Obvious candidates for such an implementation include Sine and Cosine, and, in fact, the Microsoft DirectX 9 specification does give an example using such.
One problem with only having a single function is that it may be difficult to maintain accuracy through the entire domain of the function without having a very large number of terms. A simple remedy is to again use a table approach and use a piecewise approximation to the function. This method is well known in the art and several improvements have been suggested: e.g by Detrey and de Dinechin in “Table-based polynomials for fast hardware function evaluation” or Walters and Schulte in “Efficient Function Approximation Using Truncated Multipliers and Squarers”. In this method the domain of the function is broken into 2N sections with a table entry for each section. Each entry, i, would store the necessary power series constants, i.e., C0i, C1i . . . Cki which describe the curve in that section. In practice, the approximating function would generally be linear, quadratic, or cubic. A simple example of such a quadratic scheme, with N=3, i.e. only 8 sections, is shown graphed in FIG. 1 (A). A possible look-up table for this case is shown below:
TableEntryC0C1C201.0000−0.12420.013110.8889−0.09820.009420.8000−0.07970.006930.7273−0.06590.005340.6667−0.05540.004150.6154−0.04720.003360.5714−0.04070.002670.5333−0.03550.0022Note:In a hardware implementation, if so desired, the C0i values could all be adjusted by 0.5 to reduce the cost of evaluation since this ‘bit’ is constant.
A plot of the relative error for these values is shown in FIG. 1 (B).
Given the input function value, x′, x′ε[0,1), the required table entry is identified by calculating └x′*2N┘. In hardware or software terms, this is simply a matter of taking the top N bits of a fixed-point representation of x′. An interpolating parameter, αε[0,1), is computed as α=x′*2N−└x′*2N┘. This is just a matter of extracting the remaining, least significant bits of the x′ value. The function is then computed as:ƒ(x′)=C0+C1α+C2α2 (or via the alternative factorisation presented above). This example would not be terribly accurate but can be improved through some combination of the following:
1. Increasing the number of table entries, i.e. increasing N.
2. Using a higher order polynomial.
3. Increasing the precision of the stored values.
The actual size of the multipliers and adders used in the above evaluation is also important and also represents a trade-off of implementation cost and accuracy.
Another advantage of this scheme is that, in a hardware solution, several different functions could share the same polynomial evaluation hardware, i.e., ƒ(x)=C0+C1α+C2α2+ . . . can be reused. Each particular function would, of course, still need its own look-up table and set-up and finalise units.
Although N, here, is very much smaller than that needed for the simple look-up table described earlier, it still requires a total of k2N data values. One other concern is that, because the sections are piecewise, care must be taken with the precision so that entire function is continuous, i.e., that the computed end point of one section evaluates to the matching start of the next section. It is also desirable that there not be too great a discontinuity in the derivatives across the segments.
In computer graphics, smooth surfaces and curves are frequently modelled using parametric splines or cubic Catmull-Rom. A good introduction to these systems can be found in ‘Advanced Animation and Rendering Techniques’ by Watt and Watt [ISBN0-201-54412-1]. Spline systems use ‘control points’ to control the local path of a curve or surface, which is a piecewise polynomial of some order, typically cubic. Although some splines are interpolating, i.e., the curve actually passes through the control points, in two of the more popular representation systems, Bezier and (uniform) B-Spline curves, the curve is approximating, meaning that the curve does not generally pass through the control points although it does tend to follow the positions of the control points.
In FIG. 2, a sequence of control points, (one of which is indicated by ‘1’) has been shown (connected by lines) which, in turn, defines the path of a uniform B-Spline curve, ‘2’. One particularly nice property of uniform B-Splines is that it is trivial to obtain continuity in both the curve and its derivatives across the piecewise sections. Evaluation, on the other hand, of the splines may be better performed using the de Casteljau algorithm on the Bezier equivalents. (See ‘Advanced Animation and Rendering Techniques’). Although, in graphics, splines are generally for 2D or 3D curves, the range can be reduced to a single dimension.